Load packages

library(amerifluxr) # Package for querying, downloading, and handling AmeriFlux data and metadata.
library(REddyProc)  # Package for processing Eddy Covariance flux data (e.g., gap-filling, partitioning)
library(lubridate)  # Provides tools to work with date and time (e.g., parsing, formatting).

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(tidyverse)  # A collection of packages for data manipulation, visualization, and more.

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr   1.1.4     ✔ readr   2.1.5
## ✔ forcats 1.0.0     ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1     ✔ tibble  3.2.1
## ✔ purrr   1.0.2     ✔ tidyr   1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(bigleaf)  # Provides functions to calculate leaf area index and other plant trait estimations.
library(data.table)  # Offers fast and memory-efficient manipulation of large datasets.

## 
## Attaching package: 'data.table'
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year

library(ggplot2)

01: Download AmeriFlux BASE publish

BASE Publish: After passing the Data QA/QC assessment, AMP formats the submitted flux/met data in the FP Standard format, bundles them with BADM*, and versions and publishes the resulting BASE-BADM data with a DOI (https://ameriflux.lbl.gov/data/flux-data-products/base-publish/)
AmeriFlux BASE data pipeline: https://doi.org/10.1038/s41597-023-02531-2
R package amerifluxr: https://github.com/chuhousen/amerifluxr

# check out data coverage for different sites
full.list = amf_data_coverage(data_product = "BASE-BADM", data_policy = "CCBY4.0")
head(full.list)

##   SITE_ID                                             URL
## 1  AR-Bal https://ameriflux.lbl.gov/sites/siteinfo/AR-Bal
## 2  AR-CCa https://ameriflux.lbl.gov/sites/siteinfo/AR-CCa
## 3  AR-CCg https://ameriflux.lbl.gov/sites/siteinfo/AR-CCg
## 4  AR-Cel https://ameriflux.lbl.gov/sites/siteinfo/AR-Cel
## 5  AR-Rom https://ameriflux.lbl.gov/sites/siteinfo/AR-Rom
## 6  AR-TF1 https://ameriflux.lbl.gov/sites/siteinfo/AR-TF1
##                                          publish_years
## 1                                           2012, 2013
## 2 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020
## 3             2018, 2019, 2020, 2021, 2022, 2023, 2024
## 4                                                     
## 5                                                     
## 6                                     2016, 2017, 2018

print(paste("number of sites included in BASE publish is", nrow(full.list)))

## [1] "number of sites included in BASE publish is 759"

# download single site flux/met data
site_name = "US-Syv"

# please keep the codes below commented if you do not have an account 
# floc<- amf_download_base(
#     user_id = "my_user", # change to your user name
#     user_email = "my_email@gmail.com", # change to your email
#     site_id = site_name,
#     data_product = "BASE-BADM",
#     data_policy = "CCBY4.0",
#     agree_policy = TRUE,
#     intended_use = "other",
#     intended_use_text = "Flux-LSM_workshop",
#     verbose = TRUE,
#     out_dir = tmp_dir
#   )

# If you do not have an account, you can directly start from here
tmp_dir <- tempdir()
setwd("~/") # binder wd
setwd("C:/Users/yl763/Documents/GitHub/FCC_workshop_flux_test") # Yujie's local wd, will be removed later!
floc =  "AMF_US-Syv_BASE-BADM_29-5.zip"
base <- amf_read_base(file = floc,
                      unzip = TRUE,
                      parse_timestamp = TRUE) 

# check data coverage for meteorology
# visualizes the BASE data availability for selected AmeriFlux sites, variables, and years.
##PPFD_IN   Photosynthetic photon flux density, incoming    µmolPhoton m-2 s-1
## TA   Air temperature deg C
## RH   Relative humidity, range 0-100  %
## VPD  Vapor Pressure Deficit  hPa
# more variable names: https://ameriflux.lbl.gov/data/aboutdata/data-variables/#base
amf_plot_datayear(site_set = site_name, 
                  var_set = c("PPFD_IN", "TA", "RH", "VPD"),
                  # year_set = c(2004:2023),
                  nonfilled_only = TRUE)

# check data coverage for flux data
amf_plot_datayear(site_set = site_name, 
                  var_set = c("FC", "H", "LE"),
                  # year_set = c(2004:2023),
                  nonfilled_only = TRUE)

# You can uncomment the line below to check data availability by year for all variables
# amf_plot_datayear(site_set = site_name, nonfilled_only = FALSE)

02: Organise input data

# subset data for two years
base = base[base$YEAR %in% c(2020,2021), ]

# Recreate other time variables using TIMESTAMP (This is just because I sometimes get error messages in the next step "Intialize EProc" when using the default time variables.)
recreate_time_vars <- function(df) {
  df %>%
    mutate(
      TIMESTAMP = substr(TIMESTAMP_START, 1, 12),
      year = substr(TIMESTAMP_START, 1, 4),
      month = substr(TIMESTAMP_START, 5, 6),
      day = substr(TIMESTAMP_START, 7, 8),
      hour = as.numeric(substr(TIMESTAMP_START, 9, 10)) + c(0.5, 1),
      date = as.Date(paste(year, month, day, sep = "-")),
      doy = yday(date)
    )
}
base <- recreate_time_vars(base)
# format the data to be used as input for REddyProc
base_df = data.frame(
  'Year' = base$YEAR,
  'Hour' = base$hour,
  'Date' = base$date,
  'Month' = base$month,
  'DoY' = base$DOY,
  'NEE' = base$FC,
  'LE' = base$LE,
  'H' = base$LE,
  'Rg' = ifelse(base$PPFD_IN_PI_F_1_1_1 < 0, 0, PPFD.to.Rg(base$PPFD_IN_PI_F_1_1_1)),
  'Tair' = base$TA_1_1_1,
  'Tsoil' = base$TS_1_1_1,
  'RH' = ifelse(base$RH_1_1_1 > 100, 100, base$RH_1_1_1),
  'VPD' = base$VPD_PI_F_1_1_1,
  'Ustar' = base$USTAR_1_1_1,
  'TIMESTAMP_START' = as.character(base$TIMESTAMP_START),
  'TIMESTAMP_END' = as.character(base$TIMESTAMP_END),
  'PPFD' = base$PPFD_IN_PI_F_1_1_1
)

head(base_df)

##   Year Hour       Date Month DoY       NEE        LE         H Rg      Tair
## 1 2020  0.5 2020-01-01    01   1 0.6667798 -4.066935 -4.066935  0 -12.74600
## 2 2020  1.0 2020-01-01    01   1        NA        NA        NA  0 -13.31267
## 3 2020  1.5 2020-01-01    01   1        NA -2.449585 -2.449585  0 -13.09900
## 4 2020  2.0 2020-01-01    01   1        NA        NA        NA  0 -12.97467
## 5 2020  2.5 2020-01-01    01   1 0.5290407 -2.824548 -2.824548  0 -13.11433
## 6 2020  3.0 2020-01-01    01   1        NA        NA        NA  0 -13.14767
##       Tsoil       RH       VPD     Ustar TIMESTAMP_START TIMESTAMP_END PPFD
## 1 0.6363515 87.22335 0.2930273 0.2653910    202001010000  202001010030    0
## 2 0.6407161 88.82667 0.2447449 0.2953020    202001010030  202001010100    0
## 3 0.6424356 89.99333 0.2230272 0.3732584    202001010100  202001010130    0
## 4 0.6455239 90.04335 0.2241620 0.6959635    202001010130  202001010200    0
## 5 0.6453573 89.80999 0.2268308 0.8161998    202001010200  202001010230    0
## 6 0.6477322 89.56665 0.2316199        NA    202001010230  202001010300    0

03: Intialize EProc

?filterLongRuns

## starting httpd help server ... done

#filterLongRuns : Longer runs, i.e. sequences of numerically identical values, in a series of measurements hint to problems during a noisy measurement, e.g. by sensor malfunction due to freezing. This function, replaces such values in such runs with NA to indicate missing values.
EddyData <- filterLongRuns(base_df, "NEE")
EddyData$Year <- as.numeric(EddyData$Year)
EddyData$Hour <- as.numeric(EddyData$Hour)
EddyData$DoY <- as.numeric(EddyData$DoY)
EddyDataWithPosix <- fConvertTimeToPosix(EddyData, 'YDH', Year = 'Year', Day = 'DoY', Hour  = 'Hour')

## Converted time format 'YDH' to POSIX with column name 'DateTime'.

EProc <- sEddyProc$new(site_name, EddyDataWithPosix, c('NEE','Rg','Tair','VPD', 'Ustar', "H", "LE"))

## New sEddyProc class for site 'US-Syv'

# EProc an R object used in REddyProc with the attributes defined
class(EProc)

## [1] "sEddyProc"
## attr(,"package")
## [1] "REddyProc"

04: IQR filtering

Interquartile range (IQR), which is the difference between the 75th percentile (Q3) and the 25th percentile (Q1) of a dataset;
Marginal distribution sampling (MDS, Reichstein et al., 2005) has been adopted as a standard gap-filling method;
Standard implementation of MDS relies on photosynthetic photon flux density (PPFD) (or, alternatively, global shortwave radiation) as the main driver, and air temperature (Tair) and vapor pressure deficit (VPD) as additional factors in the look-up table.
We first calculate the IQR of the residuals between MDS gap-filled values (without u* filtering) and measured values.
Here, MDS gave us an initial “best guess” of what the true flux was for each specific half hour. The choice of 6 times the IQR is a conservative approach, meaning that only the most extreme values are filtered out (some studies have used a more aggressive threshold of 3 times the IQR). We repeated the filtering process with 6 times of IQR twice to ensure sufficient outliers were removed.
One application can be found in Liu et al. (2025): https://doi.org/10.1016/j.agrformet.2025.110438

# Gapfill meteorological variables using MDS
EProc$sMDSGapFill('Tair', FillAll = FALSE,  minNWarnRunLength = NA)

## Initialized variable 'Tair' with 2810 real gaps for gap filling.

## Limited MDS algorithm for gap filling of 'Tair' with LUT(Rg only) and MDC.

## Look up table with window size of 7 days with Rg

## ............................1568
## Mean diurnal course with window size of 0 days: .
## ............0
## Mean diurnal course with window size of 1 days: .
## ............0
## Mean diurnal course with window size of 2 days: .
## ............0
## Look up table with window size of 14 days with Rg
## ............667
## Look up table with window size of 21 days with Rg
## .....575
## Finished gap filling of 'Tair' in 1 seconds. Artificial gaps filled: 35088, real gaps filled: 2810, unfilled (long) gaps: 0.

EProc$sMDSGapFill('Rg', FillAll = FALSE,  minNWarnRunLength = NA)

## Initialized variable 'Rg' with 0 real gaps for gap filling.
## Restriced MDS algorithm for gap filling of 'Rg' with no meteo conditions and hence only MDC.
## Finished gap filling of 'Rg' in 0 seconds. Artificial gaps filled: 35088, real gaps filled: 0, unfilled (long) gaps: 0.

EProc$sMDSGapFill('VPD', FillAll = FALSE,  minNWarnRunLength = NA)

## Initialized variable 'VPD' with 4270 real gaps for gap filling.
## Limited MDS algorithm for gap filling of 'VPD' with LUT(Rg only) and MDC.
## Look up table with window size of 7 days with Rg
## ..........................................2422
## Mean diurnal course with window size of 0 days: .
## ..................0
## Mean diurnal course with window size of 1 days: .
## ..................0
## Mean diurnal course with window size of 2 days: .
## ..................0
## Look up table with window size of 14 days with Rg
## ..................701
## Look up table with window size of 21 days with Rg
## ...........673
## Look up table with window size of 28 days with Rg
## ....474
## Finished gap filling of 'VPD' in 1 seconds. Artificial gaps filled: 35088, real gaps filled: 4270, unfilled (long) gaps: 0.

# Use MDS to get the "best guess" of true flux
EProc$sMDSGapFill('NEE')

## Initialized variable 'NEE' with 14679 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29511
## Look up table with window size of 14 days with Rg VPD Tair
## .......................................................582
## Look up table with window size of 7 days with Rg
## .................................................2982
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................770
## Look up table with window size of 21 days with Rg
## ...........685
## Look up table with window size of 28 days with Rg
## .....513
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 14679, unfilled (long) gaps: 0.

# Calculate residuals and identify outliers
residual <- EProc$sTEMP$NEE_orig - EProc$sTEMP$NEE_fall
IQR <- IQR(residual, na.rm = TRUE)
outlier <- ifelse(abs(residual) > (IQR * 6), 1, 0)
EddieC <- data.frame(
    sDateTime = EProc$sTEMP$sDateTime,
    NEE_orig = EProc$sTEMP$NEE_orig,
    Ustar = EProc$sDATA$Ustar,
    NEE_fall = EProc$sTEMP$NEE_fall,
    residual = residual,
    outlier = outlier
  )
  
# Rename columns
colnames(EddieC) <- c('sDateTime', 'NEE_orig', 'Ustar', 'NEE_fall', 'residual', 'outlier')
  
# Filter out outliers
EddieC$NEE_filt <- dplyr::if_else(EddieC$outlier > 0, NA_real_, EddieC$NEE_orig)
EddieC$year <- substr(EddieC$sDateTime, 1, 4)
EddieC$doy <- strftime(EddieC$sDateTime, format = "%j")
# Plot the initial outlier detection
EddieC %>%
    arrange(as.factor(outlier)) %>%
    ggplot(aes(y = NEE_orig, x = as.numeric(doy), color = as.factor(outlier))) +
    geom_point(shape = 20, alpha = 0.4) +
    theme_minimal() +
    labs(x = 'Day of year', y = 'NEE') +
    scale_color_manual(values = c('skyblue', 'red')) +
    facet_wrap(~year) + 
    ylim(c(-50, 50)) +
    ggtitle("intial outlier detection")

## Warning: Removed 14679 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Re-run the outlier test after initial filtering
EddieC$residual2 <- EddieC$NEE_filt - EddieC$NEE_fall
EddieC$IQR2 <- IQR(EddieC$residual2, na.rm = TRUE)
EddieC$outlier2 <- ifelse(abs(EddieC$residual2) > EddieC$IQR2 * 6, 1, 0)
EddieC$NEE_filt2 <- ifelse(EddieC$outlier2 == 0, EddieC$NEE_filt, NA)
  
# Plot the re-run outlier detection
EddieC %>%
    arrange(as.factor(outlier2)) %>%
    ggplot(aes(y = NEE_orig, x = as.numeric(doy), color = as.factor(outlier2))) +
    geom_point(shape = 20) +
    theme_minimal() +
    labs(x = 'Day of year', y = 'NEE') +
    scale_color_manual(values = c('skyblue', 'red')) +
    facet_wrap(~year) + 
    ylim(c(-50, 50)) +
    ggtitle("re-run outlier detection")

## Warning: Removed 14679 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Remove outliers from the main data
EProc$sDATA$NEE <- EddieC$NEE_filt2

05: u* filtering

What is u*?

u*: Friction velocity (m/s) is a measure of turbulent mixing in the lower atmosphere near the surface;
It can be thought of as a velocity scale, a representative value for a ‘turbulent velocity’;
u* is defined using the vertical flux of horizontal momentum;
(u*)^2 = overbar(u’w’)^2 + overbar(v’w’)^2;
CO2 transport is heavily influenced by turbulent mixing;
Unfavorable conditions could be detected by inspecting the relationship of nighttime NEE vs. u*.

Why is u* filtering needed?

To filter out period with low turbulence;
At low u* values, respiration is negatively biased;
At night when winds are calm, the atmosphere often becomes stable with low turbulence and low mixing locally;
This could lead to CO2 pooling near the surface or within canopies;
EC towers may underestimate respiration fluxes as CO2 is not being transported upward efficiently to be measured.

Different u* threshold treatment options:

The u* threshold is the minimum u* above which respiration reaches a plateau (see figure under /plots). - User-specific u* threshold
- Single (fixed) u* threshold
- Annually varying u* threshold
- Seasonally varying u* threshold - More details: REddyProc u* cases vignette

set.seed(2000)
# Here, we are using nSample = 10L for demonstration, please use nSample = 1000L for real research 
uStarTh <- EProc$sEstUstarThresholdDistribution(nSample = 10L, probs = c(0.05, 0.5, 0.95))

## Warning in .estimateUStarSeason(...): sEstUstarThreshold: too few finite
## records within season (n = 622) for 7 temperature classes. Need at least n =
## 700. Returning NA for this Season.

## Warning in .estimateUStarSeason(...): sEstUstarThreshold: too few finite
## records within season (n = 633) for 7 temperature classes. Need at least n =
## 700. Returning NA for this Season.

##

## Estimated UStar distribution of:
##        uStar        5%       50%       95%
## 1 0.4541172 0.4188225 0.4564736 0.5286447 
## by using  10 bootstrap samples and controls:
##                        taClasses                    UstarClasses 
##                               7                              20 
##                           swThr            minRecordsWithinTemp 
##                              10                             100 
##          minRecordsWithinSeason            minRecordsWithinYear 
##                             160                            3000 
## isUsingOneBigSeasonOnFewRecords 
##                               1

print(uStarTh)

##    aggregationMode seasonYear  season     uStar        5%       50%       95%
## 1           single         NA    <NA> 0.4541172 0.4188225 0.4564736 0.5286447
## 2             year       2020    <NA> 0.4801723 0.4488349 0.4876796 0.6412940
## 3             year       2021    <NA> 0.4280621 0.3888102 0.4258271 0.4512722
## 4           season       2020 2020001 0.4801723 0.3901247 0.4729320 0.6163569
## 5           season       2020 2020003 0.3993504 0.3877327 0.4484138 0.4839605
## 6           season       2020 2020006 0.2939407 0.2417530 0.3016328 0.3727611
## 7           season       2020 2020009 0.3376365 0.3284783 0.4120114 0.5594615
## 8           season       2021 2020012 0.4280621 0.3950345 0.4302720 0.4641630
## 9           season       2021 2021003 0.4280621 0.3742071 0.4258271 0.4512722
## 10          season       2021 2021006 0.3316664 0.2639192 0.3323809 0.3794716
## 11          season       2021 2021009 0.3632184 0.3288625 0.3656392 0.4183202
## 12          season       2021 2021012 0.4280621 0.3888102 0.4258271 0.4512722

# Define aggregation mode
EProc$sGetUstarScenarios() # by default, annual varying u* is used

##    season     uStar       U05       U50       U95
## 1 2020001 0.4801723 0.4488349 0.4876796 0.6412940
## 2 2020003 0.4801723 0.4488349 0.4876796 0.6412940
## 3 2020006 0.4801723 0.4488349 0.4876796 0.6412940
## 4 2020009 0.4801723 0.4488349 0.4876796 0.6412940
## 5 2020012 0.4280621 0.3888102 0.4258271 0.4512722
## 6 2021003 0.4280621 0.3888102 0.4258271 0.4512722
## 7 2021006 0.4280621 0.3888102 0.4258271 0.4512722
## 8 2021009 0.4280621 0.3888102 0.4258271 0.4512722
## 9 2021012 0.4280621 0.3888102 0.4258271 0.4512722

# Here, we use the single u* threshold
uStar <- uStarTh  %>%
  dplyr::filter( aggregationMode == "single") %>%
  dplyr::select( "uStar", "5%", "50%", "95%")
uStarDf <- cbind(season=na.omit(unique(uStarTh$season)), uStar)
EProc$sSetUstarScenarios(uStarDf)
EProc$sGetUstarScenarios()

##    season     uStar        5%       50%       95%
## 1 2020001 0.4541172 0.4188225 0.4564736 0.5286447
## 2 2020003 0.4541172 0.4188225 0.4564736 0.5286447
## 3 2020006 0.4541172 0.4188225 0.4564736 0.5286447
## 4 2020009 0.4541172 0.4188225 0.4564736 0.5286447
## 5 2020012 0.4541172 0.4188225 0.4564736 0.5286447
## 6 2021003 0.4541172 0.4188225 0.4564736 0.5286447
## 7 2021006 0.4541172 0.4188225 0.4564736 0.5286447
## 8 2021009 0.4541172 0.4188225 0.4564736 0.5286447
## 9 2021012 0.4541172 0.4188225 0.4564736 0.5286447

# if you want to use seasonal or annual varying u*:
# EProc$useSeaonsalUStarThresholds()# seasonal varying u*
# EProc$useAnnualUStarThresholds() # annual varying u*

# EProc$sPlotNEEVersusUStarForSeason(format = "pdf") # keep this line commented when use binder, then check the figure under /plots on github repo

Bonus training: create a figure comparing u* threshold determined using different aggregationMode.

06: MDS: gapfil NEE

The MDS combines two gap-filling techniques: the “look-up” table (LUT) and the “mean diurnal course” (MDC).

“look-up” table (LUT)

Fluxes are expected similar if they are:

at similar environmental conditions
- Rg \(\pm\) 50 \(W m^{-2}\), Tair \(\pm\) 2.5 \(°C\), and VPD \(\pm\) 5.0 \(h Pa\)
and close in time
- increasing time window until enough observations

“mean diurnal course” (MDC)

Fluxes are expected similar if they are at the same time of the day and not too many days away;
The methods are effective for short gaps as the missing values are replaced by the average of response variables under similar weather conditions in a small-time window;

Quality flag increases with fewer variables and larger time windows

0: true measurement
1: gap-filled with good quality
\(>1\): gap-filled with lower quality

# fingerplot: inspect gaps needed to be filled
EProc$sPlotFingerprintY("NEE", Year = 2020); EProc$sPlotFingerprintY("NEE", Year = 2021)

# use MDS to gapfill flux data
EProc$sMDSGapFillUStarScens('NEE')

## Ustar filtering (u * Th_1 = 0.45411720625), marked 33% of the data as gap

## Initialized variable 'NEE' with 20978 real gaps for gap filling of all 35088 values (to estimate uncertainties).

## Full MDS algorithm for gap filling of 'NEE.Ustar_uStar_fqc_0' with LUT(Rg, VPD, Tair) and MDC.

## Look up table with window size of 7 days with Rg VPD Tair

## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29068
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## .905
## Look up table with window size of 7 days with Rg
## ...................................................3092
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 20978, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.418822541913729), marked 31% of the data as gap
## Initialized variable 'NEE' with 20411 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_5%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29158
## Look up table with window size of 14 days with Rg VPD Tair
## ...........................................................843
## Look up table with window size of 7 days with Rg
## ..................................................3064
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 20411, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.456473581875), marked 33% of the data as gap
## Initialized variable 'NEE' with 21021 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_50%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29062
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## .909
## Look up table with window size of 7 days with Rg
## ...................................................3094
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 21021, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.528644729070933), marked 37% of the data as gap
## Initialized variable 'NEE' with 22024 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_95%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................28810
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## ...1076
## Look up table with window size of 7 days with Rg
## ....................................................3178
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....530
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 22024, unfilled (long) gaps: 0.

# gaps in flux data are filled 
EProc$sPlotFingerprintY("NEE_50._f", Year = 2020); EProc$sPlotFingerprintY("NEE_50._f", Year = 2021)

07: Partitioning NEE into Reco and GPP

Reco is modeled and GPP is computed as Reco - NEE;

There are three approaches to partition FCO₂ available in REddyProc:

Night-time partitioning (NT): This approach relies on the assumption that fluxes measured at night consist purely of respiration. It uses nighttime fluxes to estimate a seasonally varying reference respiration rate and its sensitivity to temperature (Reichstein 2005).
Day-time partitioning (DT): This method uses a light response curve based on daytime data to estimate the reference respiration rate, and then applies a temperature response function to estimate Reco (Lasslop 2010).
Modified day-time partitioning: Both NT and DT assume that the same Rref is applicable during daytime as at night—a questionable assumption due to the potential for inhibition of leaf respiration in the light (Keenan 2019).

# TimeZoneHour: time zone offset from UTC (Coordinated Universal Time), without daytime saving
EProc$sSetLocationInfo(LatDeg = 46.2420, LongDeg =-89.3477, TimeZoneHour = -6)
EProc$sFillVPDFromDew() # fill longer gaps still present in VPD_f
EProc$sMRFluxPartitionUStarScens("NEE") # night-time

## Start flux partitioning for variable NEE with temperature Tair_f.

## Estimate of the temperature sensitivity E_0 from short term data: 84.51.

## Regression of reference temperature R_ref for 169 periods.

## Start flux partitioning for variable NEE with temperature Tair_f.

## Estimate of the temperature sensitivity E_0 from short term data: 84.51.

## Regression of reference temperature R_ref for 169 periods.

## Start flux partitioning for variable NEE with temperature Tair_f.

## Estimate of the temperature sensitivity E_0 from short term data: 84.51.

## Regression of reference temperature R_ref for 169 periods.

## Start flux partitioning for variable NEE with temperature Tair_f.

## Estimate of the temperature sensitivity E_0 from short term data: 84.51.

## Regression of reference temperature R_ref for 169 periods.

# If you want to expolore day-time or night-time partitioning
# EProc$sGLFluxPartitionUStarScens() # day-time
# EProc$sTKFluxPartitionUStarScens() # modified day-time
names(EProc$sTEMP)

##   [1] "sDateTime"         "Tair_orig"         "Tair_f"           
##   [4] "Tair_fqc"          "Tair_fall"         "Tair_fall_qc"     
##   [7] "Tair_fnum"         "Tair_fsd"          "Tair_fmeth"       
##  [10] "Tair_fwin"         "Rg_orig"           "Rg_f"             
##  [13] "Rg_fqc"            "Rg_fall"           "Rg_fall_qc"       
##  [16] "Rg_fnum"           "Rg_fsd"            "Rg_fmeth"         
##  [19] "Rg_fwin"           "VPD_orig"          "VPD_f"            
##  [22] "VPD_fqc"           "VPD_fall"          "VPD_fall_qc"      
##  [25] "VPD_fnum"          "VPD_fsd"           "VPD_fmeth"        
##  [28] "VPD_fwin"          "NEE_orig"          "NEE_f"            
##  [31] "NEE_fqc"           "NEE_fall"          "NEE_fall_qc"      
##  [34] "NEE_fnum"          "NEE_fsd"           "NEE_fmeth"        
##  [37] "NEE_fwin"          "season"            "Ustar_uStar_Thres"
##  [40] "Ustar_uStar_fqc"   "NEE_uStar_orig"    "NEE_uStar_f"      
##  [43] "NEE_uStar_fqc"     "NEE_uStar_fall"    "NEE_uStar_fall_qc"
##  [46] "NEE_uStar_fnum"    "NEE_uStar_fsd"     "NEE_uStar_fmeth"  
##  [49] "NEE_uStar_fwin"    "Ustar_5._Thres"    "Ustar_5._fqc"     
##  [52] "NEE_5._orig"       "NEE_5._f"          "NEE_5._fqc"       
##  [55] "NEE_5._fall"       "NEE_5._fall_qc"    "NEE_5._fnum"      
##  [58] "NEE_5._fsd"        "NEE_5._fmeth"      "NEE_5._fwin"      
##  [61] "Ustar_50._Thres"   "Ustar_50._fqc"     "NEE_50._orig"     
##  [64] "NEE_50._f"         "NEE_50._fqc"       "NEE_50._fall"     
##  [67] "NEE_50._fall_qc"   "NEE_50._fnum"      "NEE_50._fsd"      
##  [70] "NEE_50._fmeth"     "NEE_50._fwin"      "Ustar_95._Thres"  
##  [73] "Ustar_95._fqc"     "NEE_95%_orig"      "NEE_95%_f"        
##  [76] "NEE_95%_fqc"       "NEE_95%_fall"      "NEE_95%_fall_qc"  
##  [79] "NEE_95%_fnum"      "NEE_95%_fsd"       "NEE_95%_fmeth"    
##  [82] "NEE_95%_fwin"      "VPDfromDew"        "PotRad_5%"        
##  [85] "FP_NEEnight_5%"    "FP_Temp_5%"        "E_0_5%"           
##  [88] "R_ref_5%"          "Reco_5%"           "GPP_5%_f"         
##  [91] "GPP_5%_fqc"        "PotRad_50%"        "FP_NEEnight_50%"  
##  [94] "FP_Temp_50%"       "E_0_50%"           "R_ref_50%"        
##  [97] "Reco_50%"          "GPP_50%_f"         "GPP_50%_fqc"      
## [100] "PotRad_95%"        "FP_NEEnight_95%"   "FP_Temp_95%"      
## [103] "E_0_95%"           "R_ref_95%"         "Reco_95%"         
## [106] "GPP_95%_f"         "GPP_95%_fqc"       "PotRad_uStar"     
## [109] "FP_NEEnight_uStar" "FP_Temp_uStar"     "E_0_uStar"        
## [112] "R_ref_uStar"       "Reco_uStar"        "GPP_uStar_f"      
## [115] "GPP_uStar_fqc"

Bonus training: comparing different partitioning method;

08: Visualise the output

FilledEddyData <- EProc$sExportResults()
combined.df <- cbind(EddyData, FilledEddyData)
names(combined.df)

##   [1] "Year"              "Hour"              "Date"             
##   [4] "Month"             "DoY"               "NEE"              
##   [7] "LE"                "H"                 "Rg"               
##  [10] "Tair"              "Tsoil"             "RH"               
##  [13] "VPD"               "Ustar"             "TIMESTAMP_START"  
##  [16] "TIMESTAMP_END"     "PPFD"              "Tair_orig"        
##  [19] "Tair_f"            "Tair_fqc"          "Tair_fall"        
##  [22] "Tair_fall_qc"      "Tair_fnum"         "Tair_fsd"         
##  [25] "Tair_fmeth"        "Tair_fwin"         "Rg_orig"          
##  [28] "Rg_f"              "Rg_fqc"            "Rg_fall"          
##  [31] "Rg_fall_qc"        "Rg_fnum"           "Rg_fsd"           
##  [34] "Rg_fmeth"          "Rg_fwin"           "VPD_orig"         
##  [37] "VPD_f"             "VPD_fqc"           "VPD_fall"         
##  [40] "VPD_fall_qc"       "VPD_fnum"          "VPD_fsd"          
##  [43] "VPD_fmeth"         "VPD_fwin"          "NEE_orig"         
##  [46] "NEE_f"             "NEE_fqc"           "NEE_fall"         
##  [49] "NEE_fall_qc"       "NEE_fnum"          "NEE_fsd"          
##  [52] "NEE_fmeth"         "NEE_fwin"          "season"           
##  [55] "Ustar_uStar_Thres" "Ustar_uStar_fqc"   "NEE_uStar_orig"   
##  [58] "NEE_uStar_f"       "NEE_uStar_fqc"     "NEE_uStar_fall"   
##  [61] "NEE_uStar_fall_qc" "NEE_uStar_fnum"    "NEE_uStar_fsd"    
##  [64] "NEE_uStar_fmeth"   "NEE_uStar_fwin"    "Ustar_5._Thres"   
##  [67] "Ustar_5._fqc"      "NEE_5._orig"       "NEE_5._f"         
##  [70] "NEE_5._fqc"        "NEE_5._fall"       "NEE_5._fall_qc"   
##  [73] "NEE_5._fnum"       "NEE_5._fsd"        "NEE_5._fmeth"     
##  [76] "NEE_5._fwin"       "Ustar_50._Thres"   "Ustar_50._fqc"    
##  [79] "NEE_50._orig"      "NEE_50._f"         "NEE_50._fqc"      
##  [82] "NEE_50._fall"      "NEE_50._fall_qc"   "NEE_50._fnum"     
##  [85] "NEE_50._fsd"       "NEE_50._fmeth"     "NEE_50._fwin"     
##  [88] "Ustar_95._Thres"   "Ustar_95._fqc"     "NEE_95%_orig"     
##  [91] "NEE_95%_f"         "NEE_95%_fqc"       "NEE_95%_fall"     
##  [94] "NEE_95%_fall_qc"   "NEE_95%_fnum"      "NEE_95%_fsd"      
##  [97] "NEE_95%_fmeth"     "NEE_95%_fwin"      "VPDfromDew"       
## [100] "PotRad_5%"         "FP_NEEnight_5%"    "FP_Temp_5%"       
## [103] "E_0_5%"            "R_ref_5%"          "Reco_5%"          
## [106] "GPP_5%_f"          "GPP_5%_fqc"        "PotRad_50%"       
## [109] "FP_NEEnight_50%"   "FP_Temp_50%"       "E_0_50%"          
## [112] "R_ref_50%"         "Reco_50%"          "GPP_50%_f"        
## [115] "GPP_50%_fqc"       "PotRad_95%"        "FP_NEEnight_95%"  
## [118] "FP_Temp_95%"       "E_0_95%"           "R_ref_95%"        
## [121] "Reco_95%"          "GPP_95%_f"         "GPP_95%_fqc"      
## [124] "PotRad_uStar"      "FP_NEEnight_uStar" "FP_Temp_uStar"    
## [127] "E_0_uStar"         "R_ref_uStar"       "Reco_uStar"       
## [130] "GPP_uStar_f"       "GPP_uStar_fqc"

# Have a look at the gapfilled data
# diurnal pattern of FC
ggplot(combined.df, aes(x = Hour, y = NEE_50._f)) +
  geom_point(col = "grey") +  
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  stat_summary(fun = mean, geom = "line") + 
  facet_wrap(~ interaction(Month, Year)) +
  theme_minimal() +
  ggtitle("FC") +  ylab(expression(FCO[2]~'('*gC~m^{-2}~y^{-1}*')'))

# diurnal pattern of GPP
ggplot(combined.df, aes(x = Hour, y = `GPP_50%_f`)) +
  geom_point(col = "grey") +  
  stat_summary(fun = mean, geom = "line") + 
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +  
  facet_wrap(~ interaction(Month, Year)) +
  theme_minimal() +
  ggtitle("GPP") + ylab(expression(GPP~'('*gC~m^{-2}~y^{-1}*')'))

## Warning: Removed 14906 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 14906 rows containing missing values or values outside the scale range
## (`geom_point()`).

# diural pattern of Reco
ggplot(combined.df, aes(x = Hour, y = `Reco_50%`)) +
  geom_point(col = "grey") +  
  stat_summary(fun = mean, geom = "line") + 
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +  
  facet_wrap(~ interaction(Month, Year)) +
  theme_minimal() +
  ggtitle("Reco") + ylab(expression(Reco~'('*gC~m^{-2}~y^{-1}*')'))

Bonus training: gap-fill H or LE, and plot the diurnal patterns.

References

Wutzler, Thomas, et al. “Basic and extensible post-processing of eddy covariance flux data with REddyProc.” Biogeosciences 15.16 (2018): 5015–5030.
Reichstein, Markus, et al. “On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm.” Global Change Biology 11.9 (2005): 1424–1439.
Liu, Yujie, et al. “Robust filling of extra-long gaps in eddy covariance CO₂ flux measurements from a temperate deciduous forest using eXtreme Gradient Boosting.” Agricultural and Forest Meteorology 364 (2025): 110438.
Keenan, Trevor F., et al. “Widespread inhibition of daytime ecosystem respiration.” Nature Ecology & Evolution 3.3 (2019): 407–415.
Lasslop, Gitta, et al. “Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation.” Global change biology 16.1 (2010): 187-208.

REddyProc

2025-08-01